Critique of Primitive Model Electrolyte Theories 
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Approximate theories for the restricted primitive model electrolyte are compared in the light of 
Totsuji's lower bound for the energy (an improvement over Onsager's), Gillan's upper bound for the 
free energy, and thermal stability requirements. Theories based on the Debye-Hiickel (DH) approach 
and the mean spherical approximation (MSA), including extensions due to Bjerrum, Ebeling, Fisher 
and Levin, and Stell, Zhou, and Yeh (PMSA1, 2, 3) are tested. In the range T* = k B TDa/q 2 < 
10 T* ~ 0.5, all DH-based theories satisfy Totsuji's bound, while the MSA possesses a significant 
region of violation. Both DH and MSA theories violate Gillan's bound in the critical region and 
below unless ion pairing and the consequent free-ion depletion are incorporated. However, the 
PMSA theories, which recognize pairing but not depletion, fail to meet the bound. The inclusion of 
excluded- volume terms has only small effects in this respect. Finally, all the pairing theories exhibit 
negative constant-volume specific heats when T* > 2T* ~ 0.1; this is attributable to the treatment 
of the association constant. 



I. INTRODUCTION 

The liquid-gas phase transition in electrolytes is of cur- 
rent interest because of puzzling experiments and theo- 
retical efforts to understand them. For recent reviews, 
see Jl)-^ . The primary model used is the restricted prim- 
itive model (RPM) consisting of two oppositely charged, 
but otherwise identical, sets of N + = N- hard spheres 
of diameter a and charge per particle ±g, immersed in 
a medium of dielectric constant D (to represent the sol- 
vent) and volume V. We will restrict our attention to 
the RPM in d = 3 dimensions and use the reduced tem- 
perature and density 



T* = k B TDa/q 2 and p* 



a 3 p, 



(1.1) 



where p = (N + + NJ)/V = N/V; the Debye inverse 
screening length, 

k d = {Aitq 2 p/Dk B T) 1/2 , with x = n D a = y/(4xp*/T*) ; 

(1.2) 

the reduced Hclmholtz free energy density 

f(p,T) = -F N (V;T)/Vk B T; (1.3) 

and the reduced configurational energy per particle, u, 
defined via 



(Nq 2 /Da) u{p, T) = U N - Mk B T, 



(1.4) 



where F N and U N = Vk B T 2 (df/dT) denote the total 
free energy and (internal) energy. 

Recent theory [^]-[| has focussed on two approaches 
to approximating the free energy of the RPM, based on 



either Debye-Hiickel (DH) theory Q or the mean spheri- 
cal approximation (MSA) [lC|-|l^]. Many years ago Bjer- 
rum fll3| ] proposed to improve DH theory by including 
ion pairing via "chemical association." Later, Ebeling 
and Grigo H] combined ion-pairing with an MSA ex- 
pression for the ionic free energy; more recently, Levin 
and Fisher || and Stell and coworkers || explored fur- 
ther extensions of the MSA. On the other hand, Fisher 
and Levin supplemented DH theory not only with 
ion pairing and excluded- volume terms but also included 
the solvation free energy of the electrically active (+, — ) 
dipolar ion pairs. Currently, this class of DH-based the- 
ories seems to give the best, albeit semiquantitative, ac- 
count of the RPM in the critical region as judged by 
comparison with simulations performed by various au- 
thors ^,|). It may be remarked that the simulation es- 
timates for T* and p* have been changing at an alarm- 
ing rate [2(b)]. Nevertheless, the MSA-based theories 
yield approximations for T* ( > 0.073) that are signif- 
icantly higher than those based on the DH approach 
(T* < 0.056), which in fact agrees much better with the 
simulations [T* = 0.048-0.056) |,|,|Jl|. 

At a purely theoretical level, however, one cannot be 
content since, a priori, there seem no clear grounds 
for preferring the DH-based theories - - apart from 
their more direct and intuitive physical interpretation — 
rather than the more modern (and fashionable) MSA- 
based theories which — since they entail the pair corre- 
lation functions and the Ornstein-Zernike (OZ) relation 
- give the impression of being more firmly rooted in 
statistical mechanics. On the other hand, it has recently 
been shown that the DH theories yield pair correlations 
satisfying the OZ relation in a very natural way |l6| ] . Fur- 
thermore, both theories have an essentially mean-field 
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character despite which, in contrast to typical mean field 
theories for lattice systems, neither has any known Gibbs- 
Bogoliubov variational formulation or similar basis. How, 
then, might the two approaches be distinguished? 

Now Blum and his coworkers have, in various places 
[fl7| [2of , enthusiastically sung the praises of the MSA for 
the RPM, asserting that the theory "is asymptotically 
correct in the limit of high density and infinite charge" or 
"high screening parameter (Debye length going to zero)." 
Furthermore, "unlike the DH theory, it [the MSA] satis- 
fies the exact Onsager bounds for the Hclmholtz free en- 
ergy and the internal energy" 0,|8) (in the same asymp- 
totic limit) and the "internal energy of the MSA is an ex- 
act lower bound" As reported below, these claims 
cannot be sustained: however, they do suggest that one 
might usefully assess and compare the MSA and DH the- 
ories, and their various extensions, by checking their pre- 
dictions against previously developed bounds for the in- 
ternal energy and Hclmholtz free energy. That task is 
undertaken here. 

Indeed, as discussed more fully in Sec. II, several 
bounds have been established. The well-known Onsager 
lower bound for the configurational energy of the RPM 
was derived in 1939 [ pT| ; less heralded is an improvement 
due to Totsuji some forty years later p2fl . For the free 
energy, Rasaiah and Stell ||23| proved that the hard-core 
free energy provides an upper bound, while Gillan p4| ] 
developed a much stronger upper bound embodying the 
idea of (+, -) pairing into dipoles [||J|J|,|l|Q. Finally, 
we note that thermodynamic stability with respect to 
temperature requires the positivity of the specific heat 
at constant volume |2q| . 

We will focus particularly on the Totsuji and Gillan 
bounds applied in the region of the predicted gas-liquid 
phase transition and critical point. We find that DH 
theory and all its augmentations always satisfy Totsuji's 
(and Onsager's) bound provided T* < 10T* ~ 0.5. On 
the other hand, the MSA actually violates the Totsuji 
bound in a significant region of the (p, T) plane where 
coexistence is predicted, unless the theory is suitably aug- 
mented. 

In the light cast by Gillan's bound, the two approaches 
rest on a more equal footing. As already shown by 
Gillan (24), the MSA (in its usual form) fails badly for 
T* < 0.08; but the same is true for the original DH the- 
ory (even when supplemented by excluded-volume terms 
) . Only when both basic theories are augmented by 
ion-pairing contributions and by allowing for the associ- 
ated depletion of the free-ion screening do they satisfy 
the Gillan bound. As against the hard-core electrostatic 
effects, included in both DH and MSA treatments, the 
presence or absence of specific excluded-volume terms 
has small effect numerically and does not affect the sat- 
isfaction of the bound. However, the recent PMSA (or 
pairing-MSA) theories of Stell and coworkers |8) violate 
Gillan's bound apparently because they do not account 
appropriately for the free-ion depletion. 

The main lesson is the crucial importance of the clus- 



tering of ions into dipolar pairs at low temperatures. Of 
course, this has been appreciated heuristically for a long 
time 1 13 and was quantitatively demonstrated in 1983 
by Gillan |2(| in calculations for the RPM which showed 
that the vapor for T* < 0.053 consisted mainly of (+, — ), 
(2+,2-), (3+,3-), ... neutral clusters and (2+, 1-), 
(1+, 2—), (3+, 2—), and (2+, 3—) singly charged clusters, 
with relatively far fewer free monopoles, (+) and (— ). 
The present work, however, seems to be the first purely 
analytic demonstration of the thermodynamic necessity 
for including clustering, implicitly or, perhaps, explicitly, 
in approximate theories. 

The recognition of (+, — ) ion-pairing requires the spec- 
ification of the corresponding association constant, K{T). 
Ever since Bjerrum's original proposal |l3f| , this has been 
a matter of confusion and contention (see, e.g., |§Q). 
Nevertheless, in the low temperature region of princi- 
pal interest here, say T* < 0.08 ~ 1.5 T*, Bjerrum's 
cutoff form and Ebeling's more sophisticated expression 
agree to within 1. 8% or better ||,[I|,|l4| and, along with 
other cutoff forms, have identical asymptotic expansions 
m powers of T* §|7). For practical purposes, there- 
fore, K{T) might be regarded as known "exactly." At 
higher temperatures, where pairing should be (and is 
predicted to be) much weaker, it is natural to surmise 
that different treatments of association would prove in- 
consequential. However, this proves false! Indeed, for 
all the previous pairing theories fl,|5|,[7|,^| Jl3|Jl^1 we find 
that the constant- volume configurational specific heat be- 
comes negative (violating thermodynamics [p5| and sta- 
tistical mechanics) in the region T* = 0.1 to 0.5: see Sec. 
IV. The source of this serious problem is found in the pro- 
posed behavior of the association constant. Initial steps 
towards amelioration are indicated, but the issue will be 
pursued in more detail elsewhere pi} ]. 

It should be mentioned that we also examine the gener- 
alized MSA (GMSA) and variants of the MSA ther- 
modynamics derived from the (approximate) pair corre- 
lation functions by routes other than the standard en- 
ergy equation p[ : these are discussed in Sec. III. Other 
even less realistic models for electrolytes exist, including 
the one-component plasma with hard cores |29| ] and the 
corresponding "dense-point limit" [11(c)]; however, we 
address here only the RPM. 

The explicit comparisons of the DH and MSA theories 
without allowance for ion pairing are presented in Sec. 
Ill, below. In Sec. IV the theories that include descrip- 
tions of ion pairing are assessed, including the PMSA 
theories M. 
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II. BOUNDS FOR THE ENERGY AND FREE 
ENERGY 

A. Configurational Energy Bounds 

The first rigorous lower bound for the configurational 
energy of the RPM seems to be due to Onsager |pl[ . It 
is essentially a consequence of the positivity of the to- 
tal electrost atic potential energy density and, with the 
notation of (L4), yields 



t(p,T) > UQns = -1- 



(2.1) 



A more transparent derivation for a system with a neu- 
tralizing background has been presented by Rosenfcld 
and Gelbart [p0[ . Totsuji, in 1981 [p2[ , improved on On- 
sager's result by writing the energy as an integral over 
the ionic pair correlation functions and showing that the 
presence of the hard-core repulsions implies an upper 
bound on the correlation functions. He thence estab- 
lished 



i(p,T)>u To t =-0.960. 



(2.2) 



Although the improvement is by only 4.0%, it has signif- 
icant consequences. 

As remarked by Totsuji, one may usefully compare 
these bounds with the electrostatic or Madelung ener- 
gies of an ionic crystal; for the NaCl (sc) and CsCl (bcc) 
structures one has 31 



"Nad 



-0.8738 and u Cs ci - -0.8813. (2.3) 



One may reasonably suppose that the latter represents 
the best possible lower bound and so we will also invoke 
it in testing approximate theories for the RPM. 



B. Gillan's Free Energy Upper Bound 



Gillan 24 has developed a convincing, but not fully 
rigorous, upper bound on the Hclmholtz free energy of 
the RPM, which incorporates the idea of ion pairing. 
The pure hard-core free energy actually provides a rig- 
orous upper bound [ p3| , but Gillan's bound is lower ex- 
cept for extremely low densities (p* < 10~ 5 ) where the 
limiting behavior is well understood. Here we utilize 
only Gillan's bound, which is derived with the aid of the 
Gibbs-Bogoliubov inequality by employing a sequence of 
truncated reference systems. The calculation finally in- 
corporates paired (+,— ) ions or dipoles by using a ref- 
erence system of over-sized, spherically-capped cylinders 
with modified Coulomb interactions. The last step of 
Gillan's argument relies on a comparison of an approxi- 
mate analytical expression for the pressure of a system of 
such spherocylinders with computer simulation estimates 
p2| , p3 |: the approximate formula appears to provide a 
bound on the true results. A search of the more recent 



literature concerning this system (e.g. Refs. |3j|-|37j]) in- 
dicates that the original simulations have withstood the 
test of time. (However, Frenkel [Q has observed that at 
high densities and for length/diameter ratios larger than 
needed here, the simulations — and, certainly, the an- 
alytic approximation — miss an isotropic-nematic fluid 
transition that is to be expected.) We thus believe that 
Gillan's bound is valid. 

To display the bound explicitly, we write the diame- 
ter and the chosen |24|] center-to-center distance of the 
spherocylinders as a s = (1 + S)a and put 

A = (5Tr/24)pal = (5tt/24)(1 + S) 3 p*. (2.4) 

If / Id (£,T) is the ideal-gas free energy density, we then 
have [E4| 



f(p,T)<-f ld (p,T) + ±pf(p,T), 



(2.5) 



T{p, T) = 1 — 2np* - -L - f A i— 1^ - In C(p, T), 



(2.6) 



C( Pl T) = T*(l — A) {1 — cxp [S/T*(l + S)}}. (2.7) 

We will adopt S = 0.3, which Gillan found optimized the 
bound for most values of T. 



III. BASIC THEORIES FOR THE RPM: 
COMPARISON WITH BOUNDS 

A. DH and MSA without pairing 

1. DH theory 

Debye-Hiickel theory |j| (here referred to as "pure" DH 
theory, since explicit dipolar pairing is not included) is 
the oldest theory for electrolytes still in current use. The 
theory entails two approximations: first, the pair corre- 
lation functions, gijiji — Vj), are represented by naive 
Boltzmann factors — with the charge, qj, multiplied by 
the average electrostatic potential at rj when an ion of 
charge qi is fixed at r.; — ignoring higher order corre- 
lation effects; and, second, these Boltzmann factors are 
linearized, which is valid only in the limit of low den- 
sity, small charge, or high temperature. (For a modern 
discussion, see McQuarrie 0.) The thermodynamics pre- 
dicted by DH theory depends only on the single parame- 
ter, x = kdo-. The appearance of the hard core diameter, 
a, demonstrates that DH theory takes account of the elec- 
trostatic effects of the hard cores; however, the original 
or pure DH theory did not treat the excluded-volume ef- 
fects of the hard cores (and so reduced to a theory for an 
ideal gas mixture in the limit of vanishing charge, q — > 0). 
Nonetheless, excluded volume contributions may be in- 
cluded naturally by adding to the free energy a suitably 
chosen pure hard-core term see below. In the DH 

critical region, such terms have a relatively small effect. 
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2. The MSA and variants 

The other "basic" theory we consider, the mean spher- 
ical approximation 10 1, is defined by a closure of the 
Ornstein-Zernike relation in which the 3y(r) vanish in- 
side the hard core, while the direct correlation functions 
outside the hard-core exclusion zone are approximated 
by the Coulombic potentials. Waisman and Lebowitz fll| ] 
solved the MSA exactly for the RPM: that is, they deter- 
mined the correlation functions which, in principle, yield 
the thermodynamics. The electrostatic free energy again 
depends only on x — kdcl, but it and the overall free 
energy depend strongly on the theoretical route taken — 
via, in particular, the energy, pressure, or compressibil- 
ity relations. Since very different results are obtained, we 
review them briefly. The standard MSA thermodynam- 
ics almost invariably discussed in the literature employs 
the energy route; but as a result, no excluded-volume 
hard-core terms are generated. Typically this problem 
is overcome by adding in appropriate terms "by hand," 
just as for DH theory |QJj|. In light of this fact, the 
conceptual advantage sometimes claimed for the stan- 
dard MSA in comparison to DH theory (see, e.g. [8(b)]), 
namely, that the former treats the hard cores in better 
fashion, seems strictly inconsequential. Note also that 
the density-density correlation functions, G pp (r), and 
also charge-charge correlation functions, G qq (r), that sat- 
isfy the Stillinger-Lovett second-momcnt-condition follow 
from DH theory (again contrary to [8(b)]) when properly 
generalized Jlq ]. 

The pressure route to MSA thermodynamics (which 
we will denote MSpA) generates a different approxima- 
tion for the electrostatic excess free energy, along with 
the Percus-Yevick-pressure-equation hard-core free en- 
ergy. It is interesting that, like the ordinary energy-route 
MSA thermodynamics, the MSpA yields both a critical 
point and the exact DH limiting laws; early on, however, 
Waisman and Lebowitz [11(c)] dismissed it as inferior. 
By contrast, the compressibility route yields no electro- 
static contribution, but generates only the Percus-Yevick- 
compressibility-equation free energy for uncharged hard 
spheres! Finally, note that the thermodynamics of the 
generalized MSA or GMSA (which is designed so that all 
three routes to the thermodynamics agree) QJl^] is iden- 
tical to the ordinary, energy-route MSA combined with 
the Carnahan-Starlin g (C S) approximation for the pure 
hard-core free energy |39fl . 



3. Hard Cores 

Since the RPM consists of hard spheres, it is certainly 
desirable to include an account of the excluded volume 
effects in any approximate theory. As we have seen, the 
two principal approximations, DH and MSA, require the 
insertion of hard cores terms "by hand," and two other 
theories, MSpA and GMSA, entail two different hard- 



core approximations. For the sake of convenience and 
uniformity, then, we will employ the CS hard core ap- 
proximation |39| in the calculations reported here for all 
theories that recognize excluded volume effects. The cor- 
responding theories will be denoted DHCS, MSACS, and 
MSpACS, while the notation DH, MSA, and MSpA will 
be reserved for the "pure" (electrostatics only) theories. 
We have, however, checked that other approximations 
for the pure hard-core contributions yield qualitatively 
similar results. 

It is worth mentioning that although hard-core terms 
do not contribute directly to the internal energy (since 
their contribution to the energy of allowed configurations 
vanishes — as correctly reflected by the CS approxima- 
tion), they do influence the overall internal energy pic- 
ture. Specifically, for the basic theories, as we shall see, 
they affect internal energy isotherms by altering the co- 
existence curve; for the augmented, pairing theories, they 
enter by changing the degree of pairing. 



B. Assessment of Basic Theories 

1, DH Configurational Energy 

For pure DH theory (with neither pairing nor hard-core 
effects) the configurational energy assumes a particularly 
simple form, namely, 



u DH (p* 



T* 



-x/2(l 



(3.1) 



Evidently the energy of DH theory viol ates n one of the 
boun ds for any values of p and T: see QL.2|), (2.1), and 



2f). Furthermore, 



,DH 



remains above the crystal val- 
ues ( [2T3| ) as is apparent in Fig. 1. The contrary state- 
ments by Blum and coworkers 0-|l9| that u DH violates 
Onsager's bound perhaps mistake the Debye-Huckcl lim- 
iting law (DHLL) — i.e., truncation of DH theory to 
lowest order in x, which no one should take seriously for 
x > 0.3: see Fig. 1 — for the full DH theory propounded 

Strictly, the d epe ndence of u on the single param- 
eter x given in (3.1) can be correct only in single- phase 
regions of the (p, T) plane. Below the critical tempera- 
ture (as defined by the theory at hand) the energy in the 
coexistence region is always a weighted sum of the values 
in the two phases, say a and (5. In fact, if the energies 
per particle are u a and up and the densities p a = p Q (T) 
and pp = pp(T), one finds 



u(p*,T*) = 



Pa{Pp - P)Ug + Pp(p — Pa)up 
p{Pf3 ~ Pa) 



(3.2) 



so that u varies linearly with 1/p. Thus the main DH 
plot in Fig. 1 is restricted to T > T® H , and similarly for 
the other theories. However, including phase coexistence 
according to (3.2) cannot induce bound violation, since a 
weighted sum of two acceptable values also satisfies the 
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bound: see the inset in Fig. 1 where the solid curves 
depict DH isotherms for T < T° H . 




FIG. 1. The configurational energy per particle for the De- 
bye-Hiickel (DH), mean spherical approximation (MSA) and 
related theories above criticality for comparison, with lower 
bounds. For a description of the bounds and the theories, see 
the text. The inset shows isotherms for T < T c for the DH 
and DHCS theories as solid and dashed curves, respectively. 
(Here and below, CS denotes use of the Carnahan-Starling 
approximation for the excluded- volume effects.) 



Regarding the effects of hard cores, one finds that the 
only changes in DHCS theory occur in the two-phase re- 
gions below T® HCS : the energy isotherms are shifted 
from those of pure DH theory since the coexistence curve 
differs. The dashed curves in the inset to Fig. 1 show the 
rather small effects: the shifts mainly reflect the expected 
lowered densities on the liquid branch of the coexistence 
curve. Naturally, these changes cannot induce any viola- 
tion of Totsuji's bound or of the crystal limits. 



2. MSA Configurational Energy 

Now Blum and Bernard have claimed the en- 

ergy of the (pure) MSA, is "asymptotically correct." 
However, as can be seen in Fig. 1, the MSA reduced 
excess energy, namely KfJ, 



,MSAi 



(p*,T*) = - [l + x - (1 + 2.x) 1 / 2 ! fx, (3.3) 



asymptotically approaches the Onsager bound of — 1 but 
violates the Totsuji bound for x > xt ~ 1200 (as Totsuji 
noted originally p2| ). Furthermore, u MSA lies below the 
crystal values for x > sx — 125. 

In fact, even in the absence of Totsuji's result, it is hard 
to make sense of the claim |lj],[l^] that the MSA energy is 
asymptotically correct for the RPM in the limit of large 
x by virtue of its approach to Onsager's bound. Agree- 
ment with a bound is hardly proof of correctness fl4l| ! 



Furthermore, the limit x — > oo at fixed density implies 
T* ~ T/q 2 — > 0; but at low temperatures, one expects 
crystalline phases to appear for p* < p* max = V% (for fee 
sphere packing) 0] and these are not described by any of 
theories under consideration. 

It is worthwhile to interpret more explicitly the val- 
ues xt and xx, where violation by the pure MSA (no 
hard cores) occurs. On the liquid side of the coexistence 
curve, xt corresponds to violation when T* < 0.012 ~ 
(0.U)T* MSA and x x corresponds to T* < 0.035 ~ 
(0.41)T* MSA . (The first violation temperature here is 
estimated with the aid of a low-temperature asymptotic 
analysis of the pure MSA coexistence curve j42| while 
the second follows directly from a numerical evaluation.) 
The solid curves in Fig. 2 demonstrate the effects. 
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FIG. 2. Comparison of the MSA energy with bounds for 
T < T c at multiples of T* = 0.01 up to T* ~ 0.0858 (solid 
curves). The dashed curve shows the T* = 0.03 isotherm for 
the GMSA for which, presumably, violations occur only at 
much lower temperatures. 

The inclusion of hard-core terms ("by hand") in the 
pure MSA changes the liquid-side coexistence curves 
more strongly than in DH theory. Thus for the MSA 
with CS terms or, equivalently, for the GMSA, the vi- 
olations shift to much lower ratios of T/T^ MSA : this 
is clearly evidenced by the dashed coexistence isotherm 

shown in Fig. 2 for T* = 0.030 ~ (0.38)T* GMSA (with 
rp*GMSA ^ 0786 [gr^ 



3. MSpA Configurational Energy 
The energy according to the MSpA is |ll[] 

nt MSpA _ _ 



1 - (1 - VT + 2x)/a 



(3.4) 



+ 2 In (1 + 2; + Vl + 2x) -2- In 4], 

which, in the single-phase region, also depends only on 
the parameter x. As evident from Fig. 1, however, this 



5 



violates the Totsuji and Onsager bounds at it — 6.5 and 
io — 7.1, respectively. These results provide ample justi- 
fication for a disparaging evaluation of the pressure-route 
thermodynamics for the MSA. For the remainder of this 
paper, we thus omit the MSpA. 

C. DH and MSA Free Energies 

In the pure theories (in which Bjerrum ion pairing is 
not explicitly included) we find that both DH theory and 
the MSA violate Gillan's free energy upper bound. The 
entire vapor branches of both coexistence curves, as well 
as both sides of the DH critical region, are in violation. 
As shown in Fig. 3, the violations remain when hard-core 
excluded volume corrections are included. The DHCS 
and GMSA treatments exhibit very similar features, for 
the low densities of interest. Note that in Fig. 3 we fol- 
low the coexist enc e prescription for the free energy cor- 
responding to ([T^) . Note also that non- violation on one 
branch of the coexistence curve (as on the GMSA liquid 
side) is at best a qualified virtue since the construction 
of the coexistence curve depends on the free energies on 
both sides. In light of these results it is clearly imperative 
to examine theories which allow for ion pairing. 

IV. ASSESSMENT OF ION-PAIRING THEORIES 

A. Bjerrum and Beyond 

To compensate for the effects of the DH linearization 
of the electrostatic Boltzmann factor, Bjerrum jl3| pos- 
tulated association of "free" ions of (residual) density p\ 
into "bound" neutral dipolar pairs of density p 2 so that 
the overall density is 

p = Pl + 2p 2 . (4.1) 

In terms of the ideal-gas free energy density f- d (pj,T) = 
Pj[l — ln(A^pj/(j)] with mean thermal de Broglie wave- 
lengths Aj(T) and internal partition functions (j{T) 0, 
we may then write the total free energy density as [|4|,[jJ 

/" = 2/7 d (i pi ) + ]{ d {p 2 ) + f Ex ( Pl , P2 ), (4.2) 

with the excess free energy density 

f BX ( P l,P2) = I HC (P1,P2) + f Im (pi) + I DI (PUP2), 

(4.3) 

where (i) f HC denotes the pure hard-core/excluded- 
volume terms, (ii) f Ion represents the electrostatic con- 
tribution of the free ions, while (iii) f DI denotes the 
dipole-ion interaction/solvation terms iQJBj. As men- 
tioned, we take here f HC to be of Carnahan-Starling 
form j3jj with the dipoles treated as effective spheres of 
diameter a 2 = 2 1 / 3 a [|l§. 
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FIG. 3. Comparison of the free energies predicted by the 
DHCS and GMSA theories in the density-temperature plane 
with Gillan's upper bound. The bound is violated below the 
solid and dashed curves, respectively. For comparison, the as- 
sociated coexistence curves with tie-lines and critical points 
are also plotted. 

Chemical equilibrium among the + and — free ions 
and dipolar pairs is imposed via the equality p 2 = 2/ii 
of the chemical potentials. If the association constant is 
defined by K{T) = A^A^/C+C- A| = C2 (see §) and 
the reduced excess chemical potentials are 

Tlf = pf x /k B T = In 7j = ~(df Ex /d Pj ), (4.4) 

with p + — p- — ipi and 7+ = 7_ = 71, then the mass 
action law states 

= K(T; Pl ,p 2 ) = K(T)2±^. (4.5) 
P+P- 72 

The optimal expression for K(T) is a matter for debate 
[Q,H — and will be discussed further below. For reference 
purposes we adopt Ebeling's form P Jl4| , |i4| which guar- 
antees an exact representation of the RPM's electrostatic 
second virial coefficient when one uses DH theory or the 
MSA (but not the MSpA) for f Ion {pi). Note that for 
T* < 0.05 ~ T* the difference between K Eb and Bjer- 
rum's original proposal, K B \ is less than 0.01%; it rises 
to 3.0% at T* = T* MSA = 0.0858, in accord with the 
Introduction. 

Bjerrum's original theory |lj| amounts to the approx- 
imation 

DHBj: j Ex ~ f Ion ~ f DH ( Xl ) with xi = Kl a, 

(4.6) 

where n\ = Airq 2 p\j DksT represents the inverse squared 
Debye length for the free ions alone, while as usual j^] , 

f DH {x) = [In (1 + x) - x + § x 2 } /47ro 3 . (4.7) 
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Friedman and Larsen jl5| later found that the predicted 
coexistence curve was unphysical. More recently, Fisher 
and Levin [|],|]|| elucidated the peculiar "banana" shape 
of the DHBj coexistence curve (see Fig. 4 below) and 
showed it became worse when excluded-volume terms 
were added as, e.g., in DHBjCS theory. However, they 
also estimated the dipole-ion solvation term as || 

f DI = p2(aai/a 3 2 T*)Q 2 (x2), x 2 = Kl a 2 , (4.8) 

Cj 2 {x) = 3 [in (1 + x + \x 2 ) - x + \x 2 ] /x 2 w x 2 /12, 

(4.9) 

where a\ = (1.0-1.3)a is the mean dipolar size, or +/— 
ion separation, while a 2 ~ 1.1619sa represents the effec- 
tive electrostatic exclusion radius ||. (Note that all the 
results given here use a\ = a and a 2 = 1.16198a.) The 
resulting DHBjDI theories lead to sensible coexistence 
curves (see Fig. 5 below) that agree fairly well with cur- 
rent simulations [5, 2(b)]. 
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FIG. 5. Comparison of BjDICS free energies, which in- 
corporate dipole-ion solvation and Carnahan-Starling ex- 
cluded-volume terms, with the Gillan bound, as in Fig. 4. 
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FIG. 4. Pure Bjerrum pairing theories tested against 
Gillan's free-energy bound. The solid and dashed "excess con- 
tours" are labeled by the magnitudes by which the DHBj and 
MSABj reduced free energies, respectively, fall below the up- 
per bound (see text). Note the associated coexistence curves 
and the unrealistic "banana" shape of the DHBj prediction 

At an earlier stage, Ebeling and Grigo ]l4j combined 
Bjerrum pairing with the MSA by replacing f DH by 

00 



f MbA {x) = 2 + 6x + 3x z - 2(1 + 2xf' 2 /12na 6 , 

(4.10) 

with x => X\ again evaluated at p\. They also added 
excluded-volume terms. The resulting MSABj and 
MSABj CS = EGA [8(b) J4] theories yield fully accept- 
able coexistence curves || but, as mentioned, the pre- 
dicted critical temperatures are significantly too high 
[2(b),5]. 



Recently, Zhou, Yeh, and Stell (ZYS) g have extended 
Ebcling's approach by using the MSA in conjunction with 
a "reference cavity theory of association" Jig] , Their 
pairing mean- spherical approximations or PMSA theories 
may be described by 



PMSA: 

f Ex = f MSA (x) 



j CS {p)+P2(T lP ) ln( 7+7 _/ 72 ), 

(4.11) 



where x — n^a is now evaluated with the total density, 
p, and / represents the single-component Carnahan- 
Starling form, evaluated at p = pi + 2p 2 (i.e., bound pairs 
are not treated as geometrically distinc t ob jects). Note 
that p 2 is here to be determined from (4.5) once K, 71, 
and 72 are specified (see below); hence p 2 is an explicit 
algebraic function of the arguments stated in (4.11). The 
use of only the total density (in place of the free ion den- 
sity pi) results in an analytically simpler, more explicit 
formulation; but, in the light of the original DH and Bjer- 
rum arguments, it seems rather unphysical since neutral 
bound pairs cannot contribute to screening in a direct 
way. Furthermore, as we will see, this approach entails a 
significant cost in accuracy. 

The specification of the PMSA may be completed by 
first noting that ZYS also adopt Ebeling's association 
constant, K Eb (T) Then, for the activity co- 

efficients, 7+ = 7_ and 72, ZYS propose three levels of 
approximation, first: 



PMSA1: 
In 7 i = 



MSA 



/dp) T = -p MSA (T,p), 



72 = 1, 
(4.12) 
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which neglects dipole-ion contributions [cf. fl4.11| )]. Sec- 
ond, dipole-ion interactions are introduced by replacing 
the approximation 72 = 1 by 

PMSA2: 

In 72 = [2(1 + + 2.x -2-Ax-x 2 } /T*x 2 , 

« -x 2 /4T*[l + 0{x)]; (4.13) 

see [8(b)], Eq. (4.11). Finally, the dumbbell-shaped hard 
cores of a dipolar ion pair are incorporated [8(a)] by using 
the CS cavity-value contact function and incrementing 
In 72 by 

PMSA3: Aln 72 =ln[2(l-77) 3 /(2-r/)], (4.14) 

where rj = np*/6. 

PMSA3 is the preferred theory of ZYS and yields 
(T*,p*) ~ (0.0745,0.0245). PMSAl and PMSA2 give 
(0.0748,0.0250) and (0.0733,0.0229), respectively. The 
T c values are still significantly higher [8(b)] than the DH- 
based estimates, namely, T* ~ 0.052-0.057 §|@' whilc 
the simulations suggest T* ~ 0.048-0.055 [2(b), 15]. 

B. Pairing Theories vs. Gillan's Bound 

Comparison of the pairing theories with Gillan's free 
energy bound is mainly encouraging. We find that theo- 
ries that incorporate association in the Bjerrum chemical 
picture, in which the free ion density is depleted by pair- 
ing (i.e., p\= p — 2^2)7 never violate the bound. Indeed, 
even the most primitive Bjerrum theories, DHBj and 
MSABj — which include neither hard-core nor dipole- 
ion interactions — satisfy Gillan's bound for all (p*,T*) 
values tested: see Fig. 4. On the other hand, all three 
PMSA theories turn out to violate Gillan's bound in sig- 
nificant regions of the (p* , T*) plane, including nearly the 
entire vapor branches of the coexistence curves. 

As regards the MSABj and DHBj theories, the more- 
or-less vertical "excess contour lines" in Fig. 4 reveal the 
magnitude of non-violation in the density-temperature 
plane: they are loci on which Gillan's upper bound ex- 
ceeds the corresponding approximate reduced free energy 
density, —fa 3 , by the indicated amounts, ranging from 
6 x 10~ 4 up to 0.1. The associated coexistence curves are 
also shown and one may notice that the excess contours 
undergo a jump in curvature on entering the correspond- 
ing two-phase region: this res ults from the coexistence 
prescription analogous to fl3.2| ). 

Fig. 5 shows the effects of incorporating dipole-ion sol- 
vation (DI) and excludcd-volumc (CS) terms. Note that 
removing the excluded- volume terms from these BjDICS 
theories produces only slight shifts in the excess contours 
at high densities and low temperatures. 

By contrast, the solid curve in Fig. 6 marks the bound- 
ary of the region inside which the PMSA3 free energy 
violates Gillan's bound. The coexistence curve is also 
shown. (Note, however, that the coexistence prescription 



was not used here to compute the violation boundary 
within the two-phase domain.) The region of violation 
found for PMSA2 is nearly identical, while that for the 
PMSAl theory is slightly larger, extending above the cor- 
responding critical point, T^ MSA1 : see the dashed curve 
in Fig. 6. 




FIG. 6. Test of the PMSA theories against Gillan's free 
energy bound. All theories fail at low temperatures and den- 
sities: see the violation boundaries, solid for PMSA3 (the 
preferred theory) and dashed for PMSAl. The coexistence 
curve and critical point are those predicted by the PMSA3. 

In conclusion, the violations of Gillan's bound found 
previously and seen here for the PMSA theories demon- 
strate convincingly that association of oppositely charged 
ions into dipoles along with a concomitant depletion 
of free ions and their screening effects is a crucial ele- 
ment in the critical-region behavior of the RPM. Gillan's 
bound also serves to highlight interesting contrasts be- 
tween DH- and MSA-based theories: the MSA coexis- 
tence curve shifts only slightly when pairing is added 
(MSABj) yet, surprisingly, violation of Gillan's bound is 
still completely avoided; the unphysical DHBj "banana" 
coexistence curve (in Fig. 4), on the other hand, imme- 
diately points to the significance of pairing, while satis- 
faction of Gillan's bound is surprising here because the 
coexistence curve is so unconvincing. 

C. Pairing Theories vs. Energy Bounds 

Testing the pairing theories against the bounds of Tot- 
suji and Onsager yields mixed results. For a window of 
temperatures that includes the critical region, namely, 
0.015 < T* < 0.5, all the theories embodying ion associ- 
ation satisfy the energy bounds. We also find a surpris- 
ing level of agreement among the various theories as to 
the value of the critical energy per particle: see Table I. 
At low temperatures, however, some of the MSA-based 
theories violate Totsuji's bound. Moreover, at moderate 
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temperatures (T* > 0.5) aZ/ofthe pairing theories violate 
fundamental thermal stability requirements (as discussed 
in the next section); for some of the approximations, this 
is also accompanied by violation of the Totsuji and On- 
sager bounds, as explained below. 

No w th e energy for a general pairing t heor y follows 
from (4.2) via th e th ermodynamic relation (1.4) and the 
mass action law (4.5), etc., which leads to 



,Ex 



{a*T ri /p*){df 1Sx ldT*) = u Ion + u DI 



a T 



— ^rJ Ex ^T) + P ju 2 (T), 



where u 2 (T) is given by 

142 (T) =T* 2 {d\nK{T)/dT*) 



(4.15) 



(4.16) 



But this can be recognized simply as the mean energy of a 
single (+, — ) bound pair since the corresponding internal 
configurational partition function for a pair is embodied 
i n th e association constant, KJT) — see the text above 
(4.4) and L evin and Fisher m. Of course, the factor 
P2\p, T) in ( 4.15[ ) is also to be determined via the law of 
mass action (4.5). For theories of the form (4.3), one can 
further write 



(4.17) 



where the "basic" expressions for the electrostatic contri- 
bution, u Ion , are now given by the natural generalizations 
of (3d) and (3^3), namely, 



u UH {p,T) = -{ Pl /p)xx/2{l + x 1 ), 



(4.18) 



,MSA i 



Go, T) = -(pi/p) [I + Xl - (1 + 2xxY /l \ I x x . 

(4.19) 

For reference, we also quote the explicit result for u DI fol- 
lowing from the treatment of Fisher a nd L evin in lea ding 
order ||). Defining ax and a 2 as in ( fi~8| ) and |h|) [||, 
one finds 



.di 



aa\ p 2 



(na 2 )' 1 



2a\ p [3 + 3na 2 + (K,a 2 ) 2 ] 



21 ' 



(4.20) 



The corresponding expressions for the PMSA theories are 
omitted for the sake of brevity. 



TABLE I. Some critical-point parameters for various theories: T c *; u c , the reduced energy per particle; x c = (Airpl/T*) 1 ' 
the (overall) Debye parameter; and, x\ a = (47rpi c /T*) 1//2 , the screening parameter. (Note that the values quoted for x c in 
correspond here to x\ c and that the Ebeling association constant |l4| was used throughout.) 





DH 


+CS 


+Bj 


+BjCS 


+BJDI 


+BjDICS 


± c 


0.0625 


0.061 3 


0.0625 


O.O6I5 


0.057 4 


0.052s 


U c 


-0.25 


-0.241i 


-0.431s 


-0.437s 


-0.444 3 


-0.453s 


Xc 


1 


0.931s 


3.013s 


3.281i 


2.466i 


2.424 


Xlc 


1 


0.931s 


1 


0.938e 


1.122 9 


0.931s 




MSA 


+CS 


+Bj 


+BjCS 


+BjDI 


+BjDICS 


T* 

± c 


0.085s 


0.078 6 


0.085s 


0.078 7 


0.082i 


0.072 3 




-0.414 2 


-0.335s 


-0.415 7 


-0.378i 


-0.444 2 


-0.414s 


x c 


2.414 2 


1.522i 


2.7213 


2.040s 


3.072 9 


2.208s 


Xlc 


2.414 2 


1.522i 


2.414 2 


1.5319 


2.4509 


1.485 




PMSA1 


PMSA2 


PMSA3 








' c 


0.073 3 


0.074s 


0.074 5 








U c 


-0.374 


-0.426 6 


-0.426s 








Xc 


1.981 4 


2.0494 


2.032 9 









1. Low Temperatures: Violation in MSA Pairing Theories 

For T* < 0.015, evaluation of u(p, T) reveals viola- 
tions of the Totsuji bound for most of the MSA theo- 
ries. The reason turns out to be literally the same as 
for the pure MSA: in the corresponding Bj, BjCS, BjDI, 
and BjDICS theories, as well as in the PMSA1 (although 
not PMSA 2 an d 3) theory, the mass-action pairing pre- 
dicted by (4.5) becomes exponentially small as T* — > 
p9[ . As a result, all these theories revert to their ion- 
only form (i.e., MSA or MS ACS) and violations occur: 



see Fig. 1. A similar depletion of pairs occurs when 
T* -> in the DHBjDI/CS (but not DHBj/CS) theories, 
and so these theories revert to the corresponding non- 
violating DH/CS theories. These results are independent 
of whether one uses the Ebeling or Bjerrum association 
constant, or any other reasonable partition- function- like 
form, as discussed below. 
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2. Moderately Low Temperatures 

In the temperature range 0.015 < T* < 0.5, which in- 
cludes T* , all the pairing theories described in the present 
study satisfy the Totsuji bound, and hence Onsager's as 
well. Fig. 7 depicts energy isotherms for the pairing the- 
ories at T* = 0.07 . The plotted isotherms have been 
cut off for large x — Krja at the hard-core packing limit, 
Pmax = Fig. 7 also shows the location of the critical 
point of the DHBjDICS theory, which may be regarded 
as a reference point in reading Table I. The table lists 
the various critical energies and Debye parameters. As 
mentioned, there is a fair measure of agreement among 
the different pairing theories regarding the energy at crit- 
icality even though other parameters vary quite strongly. 
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FIG. 7. Plots of configurational energy isotherms for var- 
ious theories at T* = 0.07, which temperature lies above all 
DH-based estimates of T* but below all MSA-based values. 
The scalloped sections of the latter isotherms thus represent 
the two-phase regions. As regards the theories, recall that 
GMSA is equivalent to MSACS and note that the "+" im- 
plies the BjDICS extensions of the basic theories. At large 
x — K,r>a, all isotherms have been cut off at the hard-sphere 
close packing density. For reference, the critical point of the 
DHBjDICS theory (where T* = 0.052 5 ) has been marked. 



3. Violations at Moderate and High Temperatures 

Violation of the energy bounds are found again, as 
mentioned above, at higher temperatures in the range 
T* > 0.5, some 6 to 10 times greater than the estimates 
for T*. The reason for this surprising fact, however, is 
quite different from the cause at low temperatures: it 
transpires, indeed, that the form of the association con- 
stant is now crucially important. 

In fact, any theory with pairing governed by Bjerrum's 
association constant violates both the Totsuji and On- 
sager bounds when T* — > 5 — and p is large enough! Once 



noticed numerically, this behavior can be und ersto od an- 
alytically by evaluating the factor u 2 (T) in ( 4.15 ) using 
Q4.16D with K = K B i(T). To that end recall, first, the 
well known fact |§ that K B ° (T) vanishes linearly, say 
as Cbj-(1 — 2T*), when T* — > i — (and remains identi- 
cally zero for T* > |). Consequently, u 2 (T) diverges to 
—00 like — t; /{I — 2T*) in this limit. However, the factor 
p?.(p, T) in (4.15) must be evaluated via the mass action 
law (F5) and is proportional to K B: >(T); this gives 



P2 

—u 2 

P 4/972 



Pill rp*2 



dK 

dT* 



8a J 72 



(4.21) 



as T* — > | — , so that p 2 — > and p\ — ► p. Note that 
the 7i(p, T), defined via (F4), depend on the theory 



under consideration. One finds that 



c B ji 



11.6: 



this is large enough so that the pairing term ( [4.2l| ) by 
itself yields a violation of Onsager's bound when (in 
DHBj theory) p* > p*™ B] ~ 0.39 or (for MSABj) 
P* > Pnl AB ^ — 0.64. However, as the other terms in 



(4.15) are also negative, violations must arise at even 
lower densities. One finds numerically, in fact, that the 
violations occur at or below p* < 0.3 in all the theories 
with pairing governed by Bjerrum's association constant. 

One expects Ebeling's choice, K Eh (T), which provides 
a match to the exact RPM second virial coefficient and 
never vanishes p 14 44] in contrast to the singular 
vanishing of K B \T)aXT* > \ — to fare better. Never- 
theless, Ebeling's association constant leads to Onsager 
and Totsuji bound violations in the region T* ~ 0.7- 1.0 
— although only in those theories which explicitly allow 
for the excluded volume effects. The PMSA3 treatment, 
furthermore, falls into this same category of violation; 
however, PMSA1 and 2 do not because the excluded- 
volume terms there do not affect the degree of pairing. 

All the violations just described turn out to be symp- 
toms of a more serious weakness of both the Bjerrum 
and Ebeling association constants, as we will now demon- 
strate. 



D. Violations of Thermal Stability 

To pursue further the origins of the Totsuji and On- 
sager energy bound violations at T* > 0.5, consider 
the energy isochores shown in Fig. 8. The two den- 
sities p* = 0.03 and 0.1 have been chosen for display 
because they bracket the critical density; similar behav- 
ior is seen at higher and lower densities. For the pure 
DH and MSA theories, included in Fig. 8 for reference 
purposes, u(p,T) rises monotonically with T: this im- 
plies a positive constant-volume configurational specific 
heat, Cy nt (p, T). (Note that outside the two-phase re- 
gion these two energy isochores are identical to those for 
DHCS and GMSA, respectively.) 

Now the positivity of the total constant-volume spe- 
cific heat is a thermodynamic necessity dictated by the 
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Second Law p5| . For a classical particle system, how- 
ever, the conhgurational contribution must be separately 
nonnegative: this follows either, thermodynamically, by 
regarding the kinetic and conhgurational degrees of free- 
dom as thermally distinct systems or, from statistical 
mechanics, by expressing Cy (p, T) as a mean-square 
energy fluctuation which is necessarily positive at finite 
positive temperatures in any nontrivial system j5Q]. 
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FIG. 8. Energy isochores at p* = 0.03 and 0.1 (note shifted 
vertical scales) for the basic DH and MSA theories and for 
various pairing theories — solid lines for those based on DH, 
dashed lines for MSA-based. Ebeling's association constant is 
employed for all plots excepting the four bracketed isochores 
for p* =0.1 labeled K Bj , which use Bjerrum's expression 
(which vanishes at T* = i). The PMSA isochores are shown 
as dot-dash curves. The plots labeled and \u® h repre- 

sent complete ion pair association (pi =0, p2 = |p), while 
U^ 3 and itf' 6 are corresponding single-pair energies implied 
by the mass-action law. Except for these plots, the isochores 
have been cut off below T* = 0.03 because by then the extrap- 
olation below T* into the two-phase regions loses all signif- 
icance. [Note that, for the approximations considered here, 



'(P,T) 



>{p,T) and u MbA (p,T) 



outside the respective two-phase regions (see Sec. IIIA.3).] 

However, a quick perusal of Fig. 8 shows that all the 
pairing theory isochores — the solid and dashed curves 
representing DH- and MSA-based theories, respectively, 
and the dot-dash plots for PMSA1 and 3 — display re- 
gions where u(p, T) decreases as T increases. In other 
words, all the paring theories predict negative constant- 
volume specific heats and violate the Second Law! 

The reason is not far to seek. In the limit of complete 
pairing (i.e., p\ =0, p 2 — \p) all the appr oximate theo- 
ries under consideration predict, via (4.15), that the en- 
ergy should be simply that of independent dipolar pairs: 
this corresponds to the plots labeled 



Fig. 8 which derive from (4.16) and the Bjerrum and 
Ebeling forms for K(T). But, as evident from the fig- 
ure, both u 2 J (T) and uf b (T) exhibit pronounced max- 
ima in the interval T* = 0.12-0.13 and then fall sharply 



as T increases, dropping below u 2 \0) 



(0) 



at T* = 0.222 and 0.21g, respectively. It is this behavior 
that leads to the decreasing regions in the overall excess 
energy isochores with incomplete pairing. But such a 
variation of u 2 (T) is physically nonsensical since, clearly, 
the configurational energy £2(r) = —q 2 /Dr of a bound 
pair cannot fall below the contact value —q 2 /Da (which, 
in turn, can be achieved in equilibrium only at T — 0). 

The problem with u 2 (T) arises because the defining 
relation ( 4.16| ) does not actually yield the physically an- 
ticipated thermodynamic mean value pl[ , say (£2 (r)) k 1 
which in the Bjerrum picture of association would be 



Mr)> 



K 



An 



e 2 (r)e' f3£2 ^r 2 dr/K(T), (4.22) 



with association constant 



r R 

K(T) = 4vr / e-^'Vdr. (4.23) 

J a 

The reason for the failure is simple: the Bjerrum cutoff 
R is taken to be temperature dependent |pl| , explicitly, 
R B 3(T) = a/2T* for T* < \ |,||. In general, such 
temperature dependence leads to the difference 



Da 



u 2 (T) 



. _ AnR 2 e- a / T ' R ^ 2 dR 
{s 2 )k ^ k B T — , 



(4.24) 



Bj 
2 d 1 



and 



which is negative whenever R(T) decreases as T rises and 
which diverges when K(T) — > 0. The Ebel ing a ssociation 
constant can also be written in the form (4.23) but with 
the large-T asymptotic form R Eb (T) -aw a/12T* 4 [§, 
which is quite accurate once T* > 0.3. We must conclude 
that neither the Ebeling nor the Bjerrum association con- 
stants can be regarded as representing even an "effective" 
partition function for an isolated ion pair as is required 
by or implicitly assumed in the standard theories of as- 
sociation urn]. 

l {p,T) As suggested by Fig. 8, the unphysically large values 

of u 2 (T) lead to negative specific heats over large re- 
gions of the (p,T) plane when either K Eb {T) or K B i (T) 
is employed. Fortunately for our primary focus on the 
critical region, the violations of thermal stability are 
confined in all cases to T* > 0.12 > 2T* (and for 
the PMSA theories to T* > 0.35). At densities below 
p* = 0.01-0.02 < 0.6p* the pairing is sufficiently weak 
that the predicted Cy (p, T) always remains positive — 
although it does display an unphysical oscillation. Once 
violations arise at a given T, moreover, they persist to 
the highest densities. 

Of course, certain features are specific to the choice 
of association constant. As remarked earlier, K B i(T) 
"switches off" abruptly at T* — i where a nonphysical 
in latent heat is implied for all p > 0; above T* — i pairing 



l„,Eb 
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is lost and no violations remain. When K Eb (T) is used 
in DH- and MSA-based theories with excluded-volume 
terms, violations remain at the highest temperatures. 

What mi ght b e a c ure for these pathologies? It is 
clear from ( 4.22 )-( 4.24 ) that the unphysical be havio r of 
u-2.(T) can be avoided if one fixes the cutoff in (4.23) at, 
say R = Xa, so defining K X (T). Furthermore, for any 
fixed A > 1, the low-T behavior of K X (T) still matches 
K Eb {T) to all orders in T* |||7). In addition, the choice 
of A may be optimized by requiring that the deviation 
\(K Eb / K x ) — 1| = 8 remain less than a specified level up 
to as high a temperature as possible. Thus one finds that 
A ~ 3.4 provides 1% precision (S = 0.01) up to T* ~ 0.11. 

One can then check that none of the pairing theo- 
ries employing K X (T) with A ~ 3.4 violates the energy 
bounds or thermal stability for any realizable thermo- 
dynamic state, (p, T) . In addition, the qualitative con- 
clusions regarding the violation and nonviolation of the 
Gillan free energy bound remain unchanged. Indeed, us- 
ing K 3A (T) causes only insignificant shifts of the free 
energy excess contours from those displayed in Figs. 4 
and 5 when T* < 0.1. 

Nevertheless, merely replacing K Eh {T) by K X (T) leads 
to significant inaccuracies in the thermodynamics at 
higher temperatures. Thus a more thoughtful approach 
is essential to providing a reasonable approximate theory 
of the RPM that is valid over the full range of tempera- 
tures [and up to moderate densities excluding, of course, 
the solid phase(s)]. Such a treatment will be presented 
elsewhere 
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